package nevo.pkgs.neuro.bio;

import static java.lang.Math.*;
import java.util.*;

import nevo.core.Utilities;

/**
 * Class that quantifies features of a spike train
 * given its voltage trace and sampling rate.
 * 
 * @author mschachter
 *
 */
public class VoltageTrace
{	
	/**
	 * Sampling rate of data (Hz)
	 */
	public long samplingRate;
	
	/**
	 * Voltage trace values (V)
	 */
	public Double[] v;
	
	/**
	 * Derivative of voltage trace (V/s)
	 */
	public Double[] dv;
	
	/**
	 * Second derivative of voltage trace in (V/s^2)
	 */	
	public Double[] d2v;
	
	/**
	 * Spike train features.
	 */
	public SpikeTrainStatistics spikeStats;
	
	
	//these parameters will not necessarily be fixed
	//in the future
	private static double vthresh = -0.05;
	private static double d2vthresh_start = -2e5;
	private static double d2vthresh_end = 0;
	private static double dvthresh_start = 0.5;
	private static double dvthresh_end = 0;
	
	private static int roffset;
	private static int foffset;
	
	
	/**
	 * Construct a VoltageTrace object given the
	 * data values and sampling rate.
	 * @param vtrace The samples of voltage (V)
	 * @param srate The sampling rate (Hz)
	 */
	public VoltageTrace(Double vtrace[], long srate)
	{	
		this.v = vtrace;
		this.samplingRate = srate;
		roffset = (int) (0.001*samplingRate);
		foffset = (int) (0.015*samplingRate);
		
		analyze();
	}
	
	/**
	 * Compute the sparse matrix that represents
	 * the phase plot, which is a plot of v vs. dv/dt
	 * 
	 * @param dvmin The minimum dv/dt (V/s)
	 * @param dvmax The maximum dv/dt (V/s)
	 * @param vmin The minimum voltage (V)
	 * @param vmax The maximum voltage (V)
	 * @param dvres The resolution of dv/dt (V/s)
	 * @param vres The voltage resolution (V)
	 * 
	 * @return A matrix representing the phase plot.
	 */
	public Double[][] computePhaseMatrix(double dvmin, double dvmax,
										 double vmin,  double vmax,
										 double dvres, double vres)
	{
		int ncols = (int) ceil(abs(vmax - vmin) / vres);
		int nrows = (int) ceil(abs(dvmax - dvmin) / dvres);
		
		Double[][] phaseMatrix = Utilities.zerosMatrix(nrows, ncols);
		
		double vavg;
		int row, col;
		for (int k = 0; k < dv.length; k++) {
			vavg = (v[k] + v[k+1]) / 2;			
			col = (int) round((abs(vmin) + vavg)/vres);
		    row = (int) round((abs(dvmin) + dv[k])/dvres);
		    
		    phaseMatrix[row][col] += 1;
		}
		
		return phaseMatrix;
	}
	
	protected void analyze()
	{
		dv = new Double[v.length - 1];
		d2v = new Double[dv.length - 1];
		
		//calculate derivatives
		dv[0] = (v[1] - v[0])*samplingRate;
		for (int k = 2; k < v.length; k++) {
			dv[k-1] = (v[k] - v[k-1])*samplingRate;
			d2v[k-2] = (dv[k-1] - dv[k-2])*samplingRate;					
		}
		
		//detect spikes
		List<Spike> spikes = new ArrayList<Spike>();
		Spike spike = null;
		int k, j, foff;
		for (k = 0; k < d2v.length; k++) {
			
			if ((spike == null) && (d2v[k] < d2vthresh_start) && (v[k] > vthresh)) {
				//a spike was detected, record threshold and time-of-arrival
				spike = new Spike();				
				for (j = k-roffset; j < k; j++) {
					//look for the start of the spike
					if (dv[j] > dvthresh_start) {
						spike.time = ((double) j)/((double) samplingRate);
						spike.thresh = (v[j] + v[j-1]) / 2;
						break;
					}
				}
				foff = min(dv.length-k, foffset);
				for (j = k; j < k+foff; j++) {
					//look for spike peak
					if ((dv[j-1] > 0) && (dv[j] < 0)) {
						spike.peak = v[j-1];
					}
				}				
			} else if ((spike != null) && (d2v[k] >= d2vthresh_end)) {
				//the end of the spike was detected, record AHP
				foff = min(dv.length-k, foffset);
				for (j = k; j < (k+foff); j++) {
					if (dv[j] > dvthresh_end) {
						spike.ahp = (v[j] + v[j-1]) / 2;
						break;
					}					
				}
				spikes.add(spike);
				spike = null;
			}
		}
		spikeStats = new SpikeTrainStatistics(spikes);
	}
}
